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The Hilbert spectrum and the Energy Preserving 
Empirical Mode Decomposition 

Pushpendra Singh, Shiv Dutt Joshi, Rakesh Kumar Patney, and Kaushik Saha 


Abstract 

In this paper, we propose algorithms which preserve energy in empirical mode decomposition (EMD), generating finite n 
number of band limited Intrinsic Mode Functions (IMFs). In the first energy preserving EMD (EPEMD) algorithm, a signal is 
decomposed into linearly independent (LI), non orthogonal yet energy preserving (LINOEP) IMFs and residue (EPIMFs). It is 
shown that a vector in an inner product space can be represented as a sum of LI and non orthogonal vectors in such a way that 
Parseval’s type property is satisfied. From the set of n IMFs, through Gram-Schmidt orthogonalization method (GSOM), n! set 
of orthogonal functions can be obtained. In the second algorithm, we show that if the orthogonalization process proceeds from 
lowest frequency IMF to highest frequency IMF, then the GSOM yields functions which preserve the properties of IMFs and 
the energy of a signal. With the Hilbert transform, these IMFs yield instantaneous frequencies and amplitudes as functions of 
time that reveal the imbedded structures of a signal. The instantaneous frequencies and square of amplitudes as functions of time 
produce a time-frequency-energy distribution, referred as the Hilbert spectrum, of a signal. Simulations have been carried out for 
the analysis of various time series and real life signals to show comparison among IMFs produced by EMD, EPEMD, ensemble 
EMD and multivariate EMD algorithms. Simulation results demonstrate the power of this proposed method. 

Index Terms 

Empirical Mode Decomposition (EMD); Energy Preserving EMD (EPEMD); Multivariate EMD (MEMD); Gram-Schmidt 
Orthogonalization Method (GSOM); Linearly independent (LI), non orthogonal yet energy preserving (LINOEP) signals. 

I. Introduction 

T he empirical mode decomposition (EMD) is an adaptive signal analysis algorithm, introduced in 11] for the analysis 
of non-stationary signals as well as signals generated from nonlinear systems and has become an established method 
for the signal analysis in various applications. The EMD decomposes a given signal into a finite number of band limited 
intrinsic mode functions (IMEs) which are derived directly from the data, unlike other signal decomposition techniques (like 
Eourier, Wavelets, etc.) which use predefined fixed basis for the signal analysis. The notion of instantaneous frequency (IE) and 
amplitude, derived from the Hilbert-Huang transform (HHT) provides an insight into the time-frequency and energy features 
of the signal. The Ensemble EMD (EEMD) is a noise-assisted data analysis method, developed in JT), to overcome the time 
scale separation problem of EMD. The Multivariate EMD (MEMD), developed in a, is a generalization of the EMD for 
multichannel data analysis. The Compact EMD (CEMD) algorithm is proposed in Ha to reduce mode mixing, end effect, and 

Pushpendra Singh, Shiv Dutt Joshi and Rakesh Kumar Patney are with the DepaiTment of Electrical Engineering, Indian Institute of Technology, Delhi, 
Hauz Khas, New Delhi-110 016, INDIA e-mails: (pushpendra.singh@ee.iitd.emet.in; sdjoshi@ee.iitd.ac.in; rkpatney@ee.iitd.ac.in). 

Kaushik Saha is with Samsung R & D Institute India - Delhi, India e-mail: (kaushik.sl4@samsung.com). 


DIGITAL SIGNAL PROCESSING 


2 


detrend uncertainty present in EMD and to reduce computation complexity of EEMD as well. To restrain the end effects and 
also to remove iterative errors and noise signal, wavelet analysis is used in the sifting process of the EMD, and a new stopping 
criterion based on correlation analysis is also proposed in ifTSl . The IMEs generated by EMD are dependent on distribution of 


local extrema of signal and the type of spline used for upper and lower envelope interpolation and the traditional EMD uses 
cubic spline for upper and lower envelope interpolation. The EMD algorithm, proposed in 113 to reduce mode mixing and 
detrend uncertainty, uses nonpolynomial cubic spline interpolation to obtain upper and lower envelopes, and have shown UB 


that it improves orthogonality among IMEs. Some recent studies, on the EMD based method have been performed for noise 
elimination m and condition-based adaptive trend prediction for rotating bearings ini. 

To eliminate energy leakage among IMEs, the Orthogonal EMD (OEMD) is proposed in 0, which generates orthogonal 
IMEs from the set of IMEs through the Gram-Schmidt orthogonalization method (GSOM). IMEs generated from EMD, EEMD 
and MEMD are not exactly orthogonal and hence there is always some energy leakage among the IME components and the 
total sum of energies of IMEs is not equal to energy of signal, i.e. energy is not preserved in decomposition. 

In any signal decomposition including EMD, the energy preserving property is important for the accurate and faithful, 
analysis and processing of three dimensional time-frequency distribution of the energy. To preserve the energy of a signal 
in the decomposition, we propose two EPEMD algorithms. Eirst EPEMD algorithm directly provides the LI, non orthogonal 
yet Energy Preserving (LINOEP) IMEs and residue (EPIMEs), and we present this novel class of ‘LINOEP’ functions in a 
well-posed mathematical result. In the second EPEMD algorithm, to completely eliminate energy leakage among the IMF 
components, we propose a method to obtain ‘orthogonal’ and ‘orthogonal & uncorrelated’ IMEs. 

This paper is organized as follows: In section II we present brief review of the various variants of the EMD (i.e. EMD, 
EEMD, MEMD) and IMEs that are required in the present work. We propose the hrst EPEMD algorithm in section III. In 
section IV, the GSOM and an orthogonal EMD (OEMD) is discussed, and we propose the second EPEMD algorithm, through 
reverse order methodology, to obtain ‘orthogonal’ and ‘orthogonal & uncorrelated’ IMEs. Simulation results are presented in 
section V. Section VI presents conclusions. 


II. The Empirical mode decomposition 


The EMD can decompose a stationary or non-stationary signal into a set of hnite band-limited IMEs. The steps involved in 
EMD algorithm Q, to extract IMEs and residue from a given signal x{t), are summarized in Algorithm 1. The sifting process 
will be continued until the hnal residue is either a constant function, or a monotonic function, or a function with only one 
maximum and one minimum from which no more IMF can be derived. The decomposed signal x{t) is expressed as the sum 
of n IMF components plus the final residue: 


n 


n+1 



( 1 ) 


i = l 


i = l 


where yi(t) is the IMF and r„(f) = is final residue. First IMF contains the finest scale or the shortest-period (i.e. 

highest frequency) oscillation and last IMF contains the longest-period (i.e. lowest frequency) oscillation present in the signal. 
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Algorithm 1: Algorithm for EMD, for i = 1, • 


1. Set yi{t) = x{t). 

2. Obtain local maxima of the signal yi{t) and generate the upper envelope e„(f) by connecting the maxima with cubic spline 
interpolation. 

3. Obtain local minima of the signal yi{t) and generate the lower envelope ei{t) by connecting the minima with cubic spline 
interpolation. 

4. Obtain the mean signal m{t) = [eu{t) + e;(t)]/2. 

5. Set yi{t) = yi{t) — m{t) and determine if yi{t) is an IMF or not by checking the properties of IMF. Repeat step 2 to 5 
and end when yi{t) is an IMF, and store it. 

6. Set x{t) = x{t) — yi{t). 

7. Repeat step 1 to 6 and end when all the IMFs and residue of signal x{t) are obtained. 

dt 


The IMFs admit amplitude-frequency modulated (AM-FMI representation ifTSll fi.e. vJt) ^ aJt) cos(6Jt)). with aJt). ^ 

0 Vt) and well-behaved Hilbert transforms m. For any IMF yi{t), its Hilbert transform yi{t) is defined as convolution of yi{t) 
and l/7rf, i.e. yi{t) = ;^ dr and the Hilbert transform emphasizes the local properties of yi{t). An analytic signal Zi{t) 

can be represented by z^{t) = yi{t) + jyi{t) = where a,{t) = [t/f(f) -I- and = tan-'^[yi{t)/yi{t)] 

are instantaneous amplitude and phase of yi{t). The IF of yi{t) is defined as; u}i{t) = 

physical meaning of IF constrains that (j)i{t) must be a mono-component function of time, and the Bedrosian and Nuttall 
theorems GUI, GD impose non-overlapping spectra constraints on the pair [ai(f),cos((/)i(f))]. An analytic representation of 
Q is given by 


i=l 

and for each IMF the IF uii{t) = For each IMF, the amplitude ai{t) and IF uii{t) are functions of time, and the 

three dimensional time-frequency distribution of amplitude is Hilbert amplitude spectrum or Hilbert spectrum 

H{uj,t) or Hilbert-Huang spectrum (HHS). The marginal spectrum which is derived from Hilbert spectrum is defined as: 
h{uj) = H(ui,t)dt. The marginal spectrum offers a measure of total amplitude (or energy) contribution from each value of 
frequency. 

The IMF components obtained from EMD methods should follow the requirements of completeness, orthogonality, locality 
and adaptiveness. The IMFs obtained from EMD satisfy the requirements of completeness, which means that the sum of the 
IMFs and residue reconstruct the original signal, and approximately follow the requirement of orthogonality IT]. All IMFs must 
satisfy two basic conditions: (1) In the complete range of time series, the number of extrema (i.e. maxima and minima) and 
the number of zero crossings are equal or differ at most by one. (2) At any point of time, in the complete range of time series, 
the average of the values of upper and lower envelopes, obtained by the interpolation of local maxima and the local minima, is 
zero. The first condition ensure that IMFs are narrow band signals and the second condition is necessary to ensure that the IF 
does not have redundant fluctuations because of asymmetric waveforms. The energy of any signal x{t), defined over the time 
[0, T], is given by x'^{t)dt and energy leakage between two IMFs can be calculated by (with j, fc = 1, 2,..., n, j A k) 

^ 3 ^ ~ An overall index of orthogonality, denoted as IOt, and a partial index of orthogonality for any two 
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IMFs components, denoted as 10jk, are defined ||T1 as follows: 


n+ln+1 pT 


•L-r i- iL-r X pi 

yj{t)yk{t)dt 

1 = 1 k=l 


IOt — 


j = L k= 
kytj 


lo x^it)dt 


rn ^ ^0 yjit)yk{t)dt 


lo y'^S)dt + \l yl{t)dt 


( 3 ) 

( 4 ) 


The ideal values of energy leakage, overall as well as partial index of orthogonality are zero. The authors in m observed that 
there is almost orthogonality among IMFs. The numerical simulations in a demonstrated that the minor eiTor in orthogonality, 
considered in HI, is not always valid, and there is actually severe energy leakage when EMD is applied for the decomposition 
of time series. The IMFs are not theoretically orthogonal, and hence the value of IOt is small but not zero, and sometimes 
very severe as shown in simulation results. The percentage error in energy (Pee) is defined as; 


Pee = x 100 


( 5 ) 


where Eemd — 's'o yi(^)dt i^ sum of energies of IMFs and residue. From ([T]), Q and Q we obtain that; 


Pee = IOt x 100 


( 6 ) 


The EMD may suffer from mode mixing, aliasing and end effect artefacts f22\ . To overcome these issues of EMD, a noise- 
assisted data analysis method EEMD, which derives the true IME components as the ensemble average of all trials, from the 
signal added with a different realization of white noise of finite amplitude in each trials. The EEMD algorithm for the signal 

x{t), to obtain IMEs, can be summarized as follows; 

Algorithm 2; Algorithm for EEMD 

1. Add a realization of white noise ni{t) ~ JV{0, to the signal x{t) i.e. set Xi{t) = x{t) + ni{t). 

2. Using EMD algorithm, obtain IMEs of the white noise added signal Xi{t) and store them. Repeat step 1 and step 2 for 
i = 1, 2, • • • ,N with different realization of white noise each time. 

3. Obtain the ensemble average of corresponding IMEs and residue of the decompositions as the final result. 

Due to ensemble averaging of corresponding IMF’s and residue, added white noise cancel each other in the hnal mean, and 

the mean IMFs stay within the natural dyadic filter windows. Therefore, EEMD preserves the dyadic property of the IMEs of 
any data, and greatly reduces the probability of mode mixing. 

The EMD and EEMD are well suited for univariate signal and when applied channel-wise to multichannel signal analysis 
may suffer form nonuniformity, scale alignment and nature of IMEs ll22l . The Multivariate EMD (MEMD) is a generalization 
of the bivariate and trivariate IT] EMD and overcomes these issues. The MEMD generates multidimensional envelopes 
by mapping multivariate signal into multiple real-valued projected signals and yield multi-dimensional rotational modes via 
the corresponding multivariate IMEs. Let, the column vectors of the matrix X = [ xi X 2 • • • x„ ] of dimension p x n, 
represent an n-dimensional multivariate signal, and the row vector d®*" = [ df' ■ ■ ■ d^ ] represent a direction 

vector corresponding to the directions given by angle 0^ = [ 9^ • • • 9^„ ] on an (n — 1) sphere, for all fc = 1 

to K directions, respectively. The MEMD algorithm, which is extension of EMD for multichannel data analysis, suitable for 
decomposition of general nonlinear and non-stationary multivariate time series is summarized from a as follows. 
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Fig. 1; Three non orthogonal vectors Ci, 02,03 such that Oi _L (02 + 03), 02 -L 03 and vector x = Oi + 02 + 03 in 3-D. 


Algorithm 3: Algorithm for MEMD 

1. Set M = X. 

2. Obtain angle 9^ corresponding to the linearly normalized low-discrepancy Hammersley sequences lfT9l on [n— 1) sphere. 

3. Obtain coordinates of unit direction vector corresponding to angle calculated in step one. 

4. Take projection of input signal on direction vector i.e. p® = d® • M^. 

5. Obtain the time instants of maxima and minima of the projected signal p® . 

6 . Generate multidimensional upper U^, and lower L^, envelopes, by spline interpolation of the input signal X at the time 
instants of maxima and minima, respectively. 

7. Obtain mean of the multidimensional upper and lower envelopes i.e. = [Ufc + L/j]/2 and set Ej, = E(fe_i) + E^, 
where Eg is null matrix. Repeat the steps 2 to 7 for k = 1 to K directions and at the end set E^, = Efe/iT. 

8 . Set M = M — Efc and determine if M is a multivariate IMFs or not. Repeat step 2 to 8 and end when M is a multivariate 
IMF, and store it. 

9. Set X = X — M. Repeat step 1 to 9 and end when all the multivariate IMFs and residue of signal X are obtained. 

The stoppage criterion for multivariate IMF's is similar to that proposed in and the hrst condition of IMF is not imposed 

as extrema cannot be properly defined for multivariate signals a. The MEMD algorithm separates out common oscillatory 
modes presents within multivariate data which makes it suitable for stationary and non-stationary multichannel data analysis. 

III. The EPEMD and LINOEP Intrinsic Mode Eunctions 
The energy preserving property is important for a variety of reasons, and it is obtained by the orthogonal decomposition of 
signal in various transforms like Eourier, Wavelet, Eourier-Bessel, etc. The energy preserving property is especially important 
for the accurate and faithful analysis and processing of three dimensional time-frequency distribution of the energy of a signal. 
The EMD algorithm, inherently, neither ensures orthogonality nor does it preserve the energy of a signal in the decomposition. 
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The IMFs, generated through EMD, are not exactly orthogonal and hence energy is not preserved (i.e. signal energy is not 
equal to sum of individual component energy). The overall index of orthogonality is not zero, and hence there is, sometimes, 
very large percentage error in the energy as shown in simulation results. To completely preserve the energy of decomposition 
and to achieve zero percentage error in energy (i.e. to obtain ideal value ‘zero’ of overall index of orthogonality), we propose 
first EPEMD algorithm to generate linearly independent (LI), non orthogonal yet energy preserving (LINOEP) IMEs and 
residue as follows. 

Applying EMD to time series x{t), one can write x{t) = yi{t) + ri(f), where, nonzero and LI signals yi{t) and ri{t) are 
the first IME and residue, respectively. 

Let sii(f) = ri{t) and Si 2 (f) = yi(f) — Q!iSii(f), where Oi = )’ sii(f) and Si 2 {t) are orthogonal. 

Through addition of Sii(<) and Si 2 {t), we obtain yi{t) + ri{t) = Si 2 {t) + (1 + ai)sii(f) = Ci{t) + C 2 (f) i.e. 

x{t) = 2/1 (f) + ri{t) = ci{t) + d^it) (7) 

where IME Ci{t) = Si 2 {t) is orthogonal to residue C 2 (f) = (1 + ai)sii(f). 

At the each stage of decomposition there are two permutations of IME and residue, and hence two choices are available to 
perform the GSOM. We proceed from the residue (lower mode of oscillations) and orthogonalize IME and residue to obtain 
proper IME Ci{t) and residue C 2 (f) which preserve the properties of the IME. If we proceed from IME (higher mode of 
oscillations) to obtain orthogonal IME and residue, then higher mode of oscillations get mixed up in residue. 

Through EMD and the above procedure, C 2 (f) and subsequent residues are further decomposed into two orthogonal components 
i.e. 


( 8 ) 


C2(i) = 2/2(f) + r 2 (t) = C 2 {t) + c'^it) ; C 2 {t) 1 C 3 (f) 
C3(^) = 2/3 (^) + = C3{t) + C4(/) ; C3(f) 1 


(9) 
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= Unit) + Vnit) = Cn{t) + C„+i(t) ; C„(<) 1 C„+i(t) 


( 10 ) 


where the residue r„(i) is not zero, if it is zero then we do not orthogonalize y„(f) and By using equations Q to 

we can write x{t) as: 

n n+1 

x{t) = YjViit) ( 11 ) 

2=1 2=1 

In this decomposition, in general, Ci{t) it Cj{t) for i,j = 1, 2,..., n, but we always have the condition that Ci{t) _L Xiji+i 
Such an example of 3D vectors is shown in Figure [T] It is to be noted that this procedure can be easily extended for MEMD, 
to obtain energy preserving MEMD (EPMEMD), and we can derive energy preserving multivariate IMEs. Based on the above 
discussions on EPEMD and LINOEP IMEs, we propose following the mathematical result. 

Theorem 1: Let iF be a Hilbert space over the field of complex numbers, and let {x,xi, • • • ,x„+i} be a set of vectors 


satisfying the following conditions: 

n + 1 

(*) X, 1 2 Xj- 

j = i + l 
n+1 

{a) X = ^ Xi 

2 = 1 

Then in the representation, given in (El, the square of the norm, and hence energy is preserved, i.e. 

2 


llxli" = 


n+1 

2 




n+1 




( 12 ) 

(13) 

(14) 


Proof: We prove the result using mathematical induction. 

Base Case: Eor n = 1, theorem is true, since x = Xi +X 2 , with Xi _L X 2 , which implies inner product (xi, X 2 ) = (x 2 , xi) = 0, 
and ||x|p = <xi + X 2 ,xi + X 2 > 

= <Xi,Xi> + <Xi,X2> + <X2,Xi> + <X2,X2> 

= I|X1|P+|1X2||1 

Induction hypothesis: Suppose that result is true for some n = A: > 1, that is, if x^ _L and x = 


Ixll" = 


Ve show that 
to show that ||x||^ = 

We have x = Xi + y, where y = Since, xi _L y, therefore, ||x|p = ll^ilP + ||y|P- From induction hypothesis we 

2 


^fc+i I 


Induction step: We show that the result is also true for n = A: + 1, that is, if x^ _L Xji= 2 >i ^ ^ hsLwe 


obtain ||y|| = 

^fc+i 


Xi = X;r =2 llxilP (It’® only a matter of indexing), and hence 


Tfc+2|| ||2 i-fi 


= X;i=i llxil! + ||xfc+ 2 |p. Hence by mathematical induction theorem is true for all positive integers n, which completes 


the proof. 


Discussion: The following observations, regarding the above result, are made: 

(1) This is a signal specific decomposition, i.e., x = 

(2) This decomposition is not unique and there may be so many sets of {x, xi, • • • , x„+i}, satisfying the above stated conditions, 
and result is valid for all such sets. 
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(3) In the context of the EMD, we come across this situation and hnd signal specific decomposition, where, we like to obtain 
signal specific IMFs which are complete, linearly independent (LI) and preserve the energy of signal. 

The EPEMD algorithm for the signal x{t), to obtain EPIMFs, is summarized in Algorithm 4: 

Algorithm 4; Algorithm for EPEMD, for z = 1, 2, • • • ,n. 

1. Set xi{t) = x{t). 

2. Using EMD algorithm, obtain IMF yi{t) and residue ri{t) of signal Xi{t). 

3. Orthogonalize IMF yi{t) and residue ri{t) (as explained above to obtain Q) and obtain new orthogonal IMF Ci{t) and 
residue 

4. Set Xi{t) = and repeat steps 2 to 4 and end when all the EPIMFs of signal x{t) are obtained. 

IV. The GSOM, orthogonal EMD and orthogonal IMFs 

In this section, we propose another algorithm, through the GSOM, by exploiting the specific order in which the properties 
of IMFs and energy of the signal are preserved in the decomposition. 

A. The GSOM and Orthogonal EMD 

The GSOM is a process for orthogonalizing a set of signals in an inner product space. Let Y = {yi{t),y 2 {t), ...,yn{t)} 
be a set of n LI signals. A set of orthogonal signals S = S 2 (f),..., s„(f)} is generated from the set Y as follows (for 

/c = 1,2,... ,n): 



yi{t) 


1 

0 . 

. 0 


Si{t) 

k-1 

Sk{t) = yk{t) - 2 Ck^s^{t) 

2 = 1 

y2{t) 

= 

C21 

1 . 

. 0 


S2{t) 


yn{t) 


Cnl 

Cn2 

. 1 


Suit) 


The Cki is obtained by using inner product (sk{t), = 0, for k ¥= i, i.e. Cki = yk(t)si(t)dt/ sj^{t)dt for i = 

1,2,... ,n, and k ^ i. where T is the total observation period of the signals. By taking sum of all the n equations of s 
along with some simple algebraic manipulations, it is shown in lO that 

n n 

= ^CiSi{t) (16) 

2=1 2=1 

where Ci = is sum of column of the coefficient matrix of ( [TSl l. It can be easily shown that Cki = 1, if fc = i. 

From set Y, there are n choices for selecting hrst signal, n — 1 choices for second signal, n — 2 choices for third signal and 1 
choice for last signal, that means there are nl permutations of the set Y, and the GSOM would produce n\ orthogonal sets of 
signals from a set of n LI signals. So it looks interesting to explore whether in a particular application, one choice is better 
than other and also why. We explore this issue in the context of EMD. 

In order to ensure the exact orthogonality and to eliminate energy leakage among IMFs, orthogonal EMD (OEMD) based 
on the GSOM is developed in |0, which generates complete orthogonal IMFs (OIMFs). Let, Y = {yi{t),y 2 {t), ...,i/„(f)} be 
a set of n IMFs of the signal x{t) generated from EMD algorithms. Through GSOM, OIMFs are obtained from a set of IMFs 
arranged in order of highest frequency to lowest frequency IMF. A signal x{t) can be expressed, in terms of IMFs, as 13 

n n 

x{t) = ^ yi{t) + Tnit) = ^ CjSi(f) + Tnit) 

2=1 2=1 


( 17 ) 
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An OEMD algorithm produces the residue signal r„(i) and n orthogonal signals Si{t), as multiplication of constant Ci on the 
signal Si{t) does not affect orthogonality. 

B. The OEMD and Orthogonal IMFs 

The OIMFs obtained from OEMD have the following limitations: (1) Higher mode of oscillation is multiplied by some 
factor and is being subtracted from lower mode of oscillation in OEMD, which would result in mixing of high frequency 
components to low frequency one, as shown in simulation results. (2) Because of the mixing of high frequency components 
to low frequency one, properties of IMFs are not preserved by signal Si{t) and hence some of the instantaneous frequencies, 
derived from the Hilbert transform, becomes negative which has no physical meaning. (3) Residue signal is not orthogonal to 
any signal component and hence there is always some energy leakage. As we have shown that the GSOM can generate n\ 
orthogonal sets from a set of n LI signals and OIMFs, generated above, are one such set. We explore other orthogonal set in 
this section. 

To overcome the limitations of OIMFs and to completely stop energy leakage, we propose second EPEMD algorithm which 
generates two sets of orthogonal IMFs through the GSOM which proceeds in the reverse order (i.e. from residue to first IMF). 
To obtain the first set of IMFs, we propose to apply the GSOM starting from the residue and finally reaching to the hrst IMF 
and we refer it as reverse orthogonal IMFs (ROIMFs). The second set of IMFs are obtained by following the same order (i.e. 
the reverse order), with only difference being that all IMFs and residue are made zero mean before applying the GSOM, and 
we refer it as reverse orthogonal and uncorrelated IMFs (ROUIMFs). 

If vectors of zero mean are orthogonal or uncorrelated, then they are ‘orthogonal and uncorrelated’ (i.e. orthogonality and 
uncorrelatedness is the same) CD, and the LI, orthogonal and uncorrelated subspaces, of inner product space, are shown in 
Figure 

The ROUIMFs are obtained through the GSOM from a set of mean removed residue and IMFs which are arranged in order 
of lowest frequency to highest frequency components (i.e. with r'^{t) = [r„(f) - rn,mean\ 

and y[{t) = \]ji{t) — yi,mean], where mean of the a signal g{t) is dehned as ^ 9{t) dt)- Thus, lower mode of oscillations is 

multiplied by some coefficient and subtracted from higher mode of oscillations in the GSOM to produces IMFs that preserve 
the properties of IMFs. It is also verihed through large number of simulations that, with this order of orthogonalization, the 
properties of IMFs are preserved. Let zi{t) = Z2{t) = y'„{t), - ■ ■ , and Zn+i{t) = y[{t). A signal x{t) is decomposed in 

n + 2 orthogonal components without any energy leakage as follows: 

n+1 

x{t) = ^ Zi{t) + C (18) 

2=1 

where constant C = Yh=i Vi,mean + Tn,mean U mean of the signal x{t). By applying the GSOM on Yli=i in ( fT^ , 

we obtain 

n+1 n+1 n+1 

= Zj = Z 

2=1 2=1 2 = 1 

The energy of the signal x{t) can be easily seen to be: 

(•T fT /n + 1 

x'^{t)dt=\ + 

Jo Jo 



(20) 
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This procedure of the GSOM produces a constant signal along with n + 1 orthogonal and uncorrelated components of the 
signal (which includes residue signal as well). There are following benefits of deriving ROUIMFs as compared to other (OIMFs 
and FOUIMFs, etc.) set of IMFs; 

(1) The ROUIMFs obtained by applying GSOM on IMFs preserve the properties of IMFs, whereas most of the FOUIMF’s 
are not able to maintain the properties of IMFs. (2) The residue signal is also orthogonal and uncorrelated to all other IMFs 
and therefore, there is no energy leakage. 

This process can be easily extended to obtain the orthogonal MEMD (OMEMD) from the MEMD, and we can derive orthogonal 

multivariate IMFs. 

We use following notations: 

{x(f)} -> EMD {yi(f), ...,t/„(t),rn(f)}={IMFs,residue}. 

{x(f)} ->EPEMD^ {ci(f),...,c„(f),c„+i(f)}={EPIMFs}. 

{yi(f), 2 / 2 (i),-..,?/nW}-GSOM^ {OIMFs} . 

{yi(f),..., t/„(f), r„(f)} ^ GSOM ^ {FOIMFsj. 

{rn{t),yn{t), ... , 2/1 (f)} ^ GSOM ^ {ROIMFs}. 

{y'i(f ),... ^ GSOM ^ {FOUIMFs}. 

{r'n{t),y'n{t). ■.. , 2 /'i(t)} - GSOM K. {ROUIMFs}. 

{xi(f), . . . , Xm{t)} MEMD {yil{t), . . . , rin{t)}, . . . , {yml{t) ; ■ • ■ ; ^mn m. 

{rin{t), ...,yii{t)}, ...,{rm„(t), ...,ymi(f)} GSOM {ROIMEs},...,{ROIMEs}. 

V. Simulation results 

The online available MATLAB software for EMD and EEMD 12^ . and for MEMD ||24| have been used in simulation 
results. The objectives of the simulations are (1) to calculate and compare energy leakage, percentage energy error and index 
of orthogonality from the proposed algorithms as well as from the EMD, EEMD and MEMD algorithms, by using the simulated 
and real file signals. ( 2 ) to use the proposed algorithm for time-frequency analysis of a chirp signal and compare result with 
the EMD and EEMD. (3) to test statistical significance of IMEs generated by the proposed algorithms. 

A. The comparison of energy leakage between EMD and EPEMD 

The overall index of orthogonality (IOt) is shown in Eigurej^for the following different type of signals, with Ai=100, 
^2=1, sampling frequency 1^8=150 and time duration of 0 to 10 sec. (1) Low pass signal, LP = [A2sm(27r(50 — i)t) + 

Aism(27r(l -I- i)t)]; (2) Band pass signal, BP = 2^°^[A2Sin(27r(50 — i)t) + Aism(27r(15 + i)t) + A2sm(27r(l -I- (3) 

High pass signal, HP = 2^°^[Aism(27r(50 — i)t) + A2sm(27r(l + i)t)]; (4) Band stop signal, BS = 2]i£i[^is*n(27r(50 — 
i)t) + A2sin(27r(15 + i)t) + Aism(27r(0 + (5) All pass signal, AP = 2^°j^[Aism(27rit)]; (6) AM signal, AM = 

(1 + A2sm(27r3f)) • (Aism(27r20f)); (7) EM signal, FM = Aism((27rl0 + 5sin(2n3t))t); (8) White Gaussian noise, WGN 
(normal distribution with mean 0 and standard deviation 1); and (9) CHIRP (linear chirp of amplitude 100 and frequency 0.1 Hz 
to 50 Hz). Eigurej^and Eigurej^show IOt for different sampling rates, generated from signal s{t) = Xi/°=i[100' sin{2Trfit)] 
with time duration of 0 to 10 sec, sampling rate Fg=lQ5 to 400 Hz with increment of 5 Hz and Fs=105 to 2000 Hz with 
increment of 50 Hz, respectively. As shown in these simulation, IOt is varying with type of signal and sampling frequency 
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Fig. 3: Overall Index of Orthogonality for different type of signal 



Fig. 4: IOt for different sampling frequency 


and sometimes it is very high (e.g. in Figure |^peak value of IOt ^ —58) in EMD whereas in case of EPEMD IOt is always 
in the range of 10“^® which is almost zero. 

B. Real life time series decomposition 

In this simulation we take real life time series to compare the energy leakage and percentage energy error among the IMEs 
generated by the EMD, EEMD and OEMD algorithms. 
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Fig. 5: IOt for different sampling frequency 


TABLE I: The overall Index of Orthogonality IOt 



IMEs 

OIMEs 

EOIMEs 

ROIMEs 

EMD 

48.5e-3 

79.2e-3 

92.3e-17 

-41.9e-18 

EEMD 

-15.1e-3 

-69.7e-3 

11.5e-17 

-56.0e-18 


TABLE II; The signal energy = 11.3752, components and sum of components energies and Pee for IMEs, OIMEs, EOIMEs 
and ROIMEs obtained from EMD. E* is energy of residue component. 



IMEs 

OIMEs 

EOIMEs 

ROIMEs 

El 

707.3e-3 

632.1e-3 

1.15 

694.7e-3 

E 2 

475.2e-3 

407.2e-3 

718.53e-3 

429.6e-3 

Ei 

277.05e-3 

115.3e-3 

601.6e-3 

195.1e-3 

Ei 

783.5e-3 

739.2e-3 

1.4 

552.2e-3 

Es 

88.5e-3 

88.3e-3 

660.52e-3 

12.7e-3 

Ee 

2.5e-6 

2.1e-6 

5.0 

73.8e-3 

E* 

8.5 

8.5 

1.8 

9.4 

E'jp 

10.8 

10.5 

11.4 

11.4 

Pee 

4.85 

7.9 

124.9e-15 

46.8e-15 


1) The annual mean global surface temperature anomaly time series analysis: In order to compare and demonstrate the 
advantage of the OEMD with ROIMEs, we applied the proposed method to data “the annual mean global surface temperature 
anomaly”, online available ||23|, as shown in Eigure The sets of IMEs, EOIMEs, ROIMEs and ROUIMEs obtained by 
the decomposition of the data through the EMD and EEMD methods are shown in Eigures to 14 which shows significant 
improvements in ROIMEs and ROUIMEs as compare to EOIMEs. The overall Index of Orthogonality IOt for EMD and EEMD 
are given in Table 11] which clearly indicates better performance of ROIMEs over others. The value of partial orthogonality 
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The annual mean global surface temperature anomaly 



Fig. 6 : Annual mean global surface temperature anomaly 


IMPS 



Fig. 7; IMFs yi to j /4 and (j /5 + ye + tq) obtained from EMD. 


index of IMFs are in the range of 10“^ and for FOIMFs and ROIMFs are in the range of 10“^®. The index of total and partial 
orthogonality, and Pee are least in ROIMFs. The value of partial orthogonality index and Pee for ROUIMFs is almost the 
same as the value of ROIMFs. 

2) The Elcentro Earthquake May 18, 1940 North-South Component time series analysis: The Elcentro Earthquake data 
has been taken from ESil and is shown in Eigure The percentage errors in total signal energies {Pee) for IMEs, OIMEs, 
and ROIMEs obtained from EMD, for the Elcentro Earthquake time series data, are given in Table [III] which indicates better 
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FOIMFs 



Fig. 8 ; FOIMFs yi to 1/4 and (1/5 + ye + rg) obtained from EMD. High frequency components are mixed in low frequency 
ones. 



1850 1900 1950 2000 

year 


Fig. 9; ROIMFs yi to 1/4 and ( 1/5 + ye + r^) obtained from EMD 


performance (i.e. least value of Pee) of ROIMEs over others. The Hilbert marginal spectrum for IMEs, EOUIMFs, ROUIMEs 


and EPIMEs derived through EMD of the ElCentro earthquake time series data is shown in Figure 16 As there are no Pee 


and no energy leakage among ROUIMFs, and no Pee in EPIMFs, therefore, the Hilbert marginal spectrum of ROUIMFs and 
EPIMFs can more accurately and faithfully characterize the signal energy distribution at each frequency components. 
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1850 1900 1950 2000 

year 

Fig. 10: ROUIMFs yi to and ( 1/5 + j/e + ^e) plus DC component obtained from EMD 

IMFs 



1850 1900 1950 2000 


year 

Fig. 11: IMFs yi to 1/4 and ( 1/5 + ye + re) obtained from EEMD 

TABLE III: Pee for IMEs, OIMEs, and ROIMEs obtained from EMD for Elcentro Earthquake data, 

IMFs OIMFs ROIMFs 

Pee -26 -2.3e-2 5.55e-13 
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Fig. 14; ROUIMFs j/i to 1/4 and (j /5 + ye + r^) plus DC component obtained from EEMD 


Elcentro Earthquake May 18, 1940. North-South Component 



Eig. 15; Elcentro Earthquake May 18, 1940 North-South Component 


for the simulation results shown in Eigure to 
decomposition. 

D. The Time-Frequency Analysis of chirp signal 

Eigure [20| shows the Time-Erequency (T-E) estimates for linear chirp (time 0 to 0.3 second, Fs = 10000 Hz, frequency 
range 100 to 200 Hz, 50 zero padded), obtained using the EMD, EEMD and EPEMD. There is enhanced T-E tracking when 
using EPEMD as compare to EMD and EEMD. The reason for the artifacts in EMD and EEMD is high energies in the IMEs 


19 The ROIMEs preserve properties of IME and energy of signal in 
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Fig. 16: Hilbert marginal spectrum for Elcentro Earthquake May 18, 1940 North-South Component time series 
TABLE IV: Pee for IMEs, EOIMEs and ROIMEs obtained from MEMD 



IMFs 

FOIMFs 

ROIMFs 

Pee xi{t) 

23.7398 

1.1505e-14 

5.7526e-14 

Pee X 2 {t) 

20.2972 

-2.9996e-14 

2.9996e-14 

Pee X 2 {t) 

29.9836 

-7.5779e-14 

-4.5467e-14 

Pee X 4 (t) 

19.0329 

-4.5440e-14 

-1.5147e-14 



Eig. 17: 4-variate signal Xj{t) (first row) and its IMEs (second row onwards) obtained from MEMD of (|2T]l in each column. 
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Fig. 18; FOIMFs obtained from MEMD. High frequency components are mixed in low frequency ones. 
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Fig. 19: ROIMFs obtained from MEMD 


of low frequencies due to leakage. The percentage energy leakage (Pee) for EMD, EEMD and EPEMD are —1.014 x 10^, 
113.54 and 3.181 x 10“^^, respectively. 

E. The statistical significance of IMFs generated by the proposed EPEMD algorithms 

The statistical significance test of IMEs is developed in lITOl to determine if data or its IMP components contain relevant 
and useful or not so relevant information. The IMF components with their energy located within the upper and lower bounds 
are considered as components generated from the Gaussian white noise present in data and contain no signal information, and 
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Time (s) 
(c) EPEMD 



Time (s) 


Fig. 20; Time-Frequency Analysis of linear chirp with zero padding. 


the IMF components with their energy located outside the upper and lower bounds contain signal information, at the selected 
confidence level. 

To study the statistical characteristics of white noise using EMD, the Fourier spectra and mean periods estimation, the 
statistical significance test for each IMFs has been performed in iflOl . Similarly, here we obtain Figures |^to 25 which are 
the spread function plot, obtained via EMD and EPEMD, for statistical significance test of the IMEs, EPIMEs, EOIMEs and 
ROUIMEs of the Gaussian white noise of 2^® samples with the normal distribution of mean 0 and standard deviation 1. The 
dashed and solid lines are the 5*^ and 95*^ percentiles, respectively. The stars correspond to the pairs of the averaged mean 
energy density and the averaged mean period of IMEs. Since we have obtained these plots from the decomposition of the 
Gaussian white noise, we expect averaged mean energy density well within confidence limit. It is clear from the figures that 
the averaged mean energy density of all EPIMEs, ROIMEs and ROUIMEs are well within confidence limit and represent the 
Gaussian white noise, whereas the averaged mean energy density of two IMEs and many EOIMEs components are not within 
confidence limit and hence these components represent spurious signal components and they are not providing any physical 
meaning. 


VI. Conclusions 

In this paper, we have proposed two energy preserving EMD (EPEMD) algorithms. The first EPEMD algorithm, to preserve 
the energy of a signal in decomposition, decomposes a signal into the linearly independent (LI), non orthogonal yet energy 
preserving (LINOEP) IMEs and residue (EPIMEs). In the second algorithm, we have shown that if orthogonalization process 
proceeds from the lowest frequency component to highest frequency IME, the GSOM yields functions which preserve the 
properties of IMEs as well as the energy in decomposition, and hence eliminate the energy leakage among IMEs. The suitability 
of the generated reverse orthogonal IMEs (ROIMEs) and reverse orthogonal and uncorrelated IMEs (ROUIMEs) are validated 
through the decomposition of various simulated as well as real life time series. The overall and partial index of orthogonality 
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Significance of IMFs 



Fig. 21: Significance test of the IMFs of the Gaussian white noise. 


Significance of EPIMFs 



Fig. 22: Significance test of the EPIMFs of the Gaussian white noise. 

and energy leakage are used to demonstrate the marked improvement in the orthogonality of the ROIMFs and ROUIMFs 
components. Finally, the statistical significance test of IMFs, EPIMFs, ROIMFs and ROUIMFs generated from the white 
Gaussian noise via EMD and EPEMD algorithms are used to illustrate the relevance, improvement and effectiveness of the 
proposed methodologies. The EPIMFs, ROIMFs and ROUIMFs generated by various EMDs preserve the properties of IMFs 
and the energy of signal in the decomposition, and are promising and generating better results. 
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Significance of FOIMFs 



Fig. 23: Significance test of the FOIMFs of the Gaussian white noise. 


Significance of ROIMFs 



Fig. 24: Significance test of the ROIMFs of the Gaussian white noise. 
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